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CN Abstract 

The fast marching algorithm, and its variants, solves numerically the generalized eikonal 
equation associated to an underlying riemannian metric A4. A major challenge for these 
algorithms is the non-isotropy of the riemannian metric, which magnitude is characterized by 
ff") the anisotropy ratio n(M) € [1, oo]. Applications of the eikonal equation to image processing 

CN [U [3] often involve large anisotropy ratios, which motivated the design of new algorithms. 

A variant of the fast marching algorithm, introduced in [6j, addresses the problem of 
large anisotropics using an algebraic tool named lattice basis reduction. The numerical 
y-— complexity of this algorithm is insensitive to anisotropy, under extremely weak assumptions. 

We establish in this paper, in the simplified setting of a constant riemannian metric, that 
the accuracy of this algorithm is also extremely robust to anisotropy : in an average sense, 
it does not degrade as n(M) increases. We also extend this algorithm to higher dimension. 

a 

Introduction 

The Generalized Eikonal Equation is a Partial Differential Equation (PDE) , which characterizes 
the riemannian distance, associated to a riemannian metric Ai, between a domain's boundary 
and an arbitrary point of this domain. Alternatively this PDE is also the level set formulation of 
an elementary front propagation model [9], where the front speed is dictated locally by the front 
position and orientation, a dependence encoded in the riemannian metric M, but independent 
^vq of global properties of the front, or higher order properties such as its curvature. This restrictive 

setting allows to compute numerically the front evolution using the fast marching algorithm or 
its variants, which are derivatives of Dijkstra's shortest path algorithm, instead of more costly 
level set methods. Solutions to the Generalized Eikonal equation have numerous applications |9 , 
including medical image processing [HE] which motivates this work. In this context, pronounced 
anisotropics are not uncommon, which challenges currently available algorithms. 

The original fast marching algorithm [10J is limited to isotropic metrics. Several variants 
were later developed, which can handle limited [3] or arbitrary [8j [6] anisotropics. Anisotropy 
magnitude is characterized by the anisotropy ratio k(A4), see ([3]), which reflects the local distor- 
tion of distances. The method [6| was developed for applications based on a cartesian grid, such 
as image processing, where robustness to large anisotropics is crucial, both in terms of compu- 
tational complexity and numerical accuracy. This algorithm was originally limited to two and 
three dimensional domains, but we extend it to four dimensions in the last section of this paper. 
The complexity 0(NlnN + N \n k(M)) of this algorithm, where N denotes the cardinality of 
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the discrete computational domain, is similar to that of the original isotropic fast marching 
algorithm 0(N\nN), under the weak assumption that the number of discretization points iV 
exceeds the anisotropy ratio k(J\A) (in typical applications k(A4) < 100 and N > 10 4 ). 

We study in this paper the accuracy of this algorithm, motivated by numerical experiments, 
which showed significant improvements over an alternative method [2j for solving the generalized 
eikonal equation. More precisely, the highly anisotropic benchmark (k(M) = 100) presented on 
Figure 3 in [B] shows an error reduction by a factor 7 on a 1200 x 1200 grid, while the compu- 
tation time is reduced by a factor 200. The error analysis in the case of an arbitrary continuous 
riemannian metric M on a general domain was not succeeded, and we therefore restrict our 
attention to case of a constant metric M on the domain lR d \ {0}. While this simplified setting 
is purely academic, the author believes that it gives a valuable insight on the local behavior of 
the general case. Note that the fast marching algorithm, and most of its variants, are first order 
discretizations of the eikonal equation PDE, but that we are interested the precise dependence 
of the numerical error with respect to the anisotropy ratio k(J\A). In the worst case scenario, 
our error bound grows like a power of k(J\A). We show however in Theorems [l] and [2] that, in 
an average sense over all space directions, the numerical error in independent of the anisotropy 
ratio k(M). 

In order to state our results, we need to introduce some notations. We consider a fixed 
integer d > 1, and we denote by Si" the collection of d x d symmetric positive definite matrices. 
For each M £ Sj" we introduce the scalar product (-,-)m and the norm || • \\m defined for all 
u, v £ TR, d by 

{u,v) M ■= u T Mv, \\u\\ M ■= y/{u,u) M - (1) 

Consider an open domain SI C ]R rf , and a riemannian metric Ai £ C°(S7, Sf). The generalized 
eikonal equation is the following PDE : 

I vd ( 2: )IIa4(^)- 1 = 1 for all z £ S7, 

D(z) = for all z £ <9S7. [ ' 

We chose null boundary conditions for simplicity, and we refer the reader interested in more 
general boundary conditions to the discussion in [2j. The PDE ([2]) characterizes the distance 
function D £ C°(S7,H + ) to the boundary of SI : the unique viscosity solution [3] is 

D(z) = inf{length( 7 ); 7 € C\[0, 1],0), 7 (0) = z, 7 (1) £ 

length( 7 ) := / \\i {t)\\ M{l{t)) dt. 
Jo 

The anisotropy ratio n{M) £ [1, oo] is the supremum value of |M|.mO)/IMI.M(,z)> where z £ SI 
and u, v £ lR d are vectors of unit euclidean norm. In other words 

k{M) := sup k(M(z)), where for M £ Sj, k(M) := ^\\M\\\\M' l \\. (3) 

The discretization of the PDE ^ takes the form of a fixed point problem |10t [2l |6] 

d(z) = A(d, z) for all z £ SI*, 
d(z) = for all z £ an*, 



(4) 



where S7* and 5S7* denote discrete sets devoted to the sampling of the continuous domain S7 and 
its boundary <9S7 respectively, and d : SI* U <9S1* — > H + denotes a discrete map. The Hopf-Lax 
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Figure 1: Three M-reduced meshes, produced using different strategies, for a matrix MeSj 
of anisotropy ratio k(M) = 10 and eigenvector (1,0.6). 



update operator A(d, z) can take several forms, but depends only on the value of d(x) for a finite 
number of points of x G f2* U<9f2*, referred to as the neighbors of z. Under suitable assumptions, 
the fixed point problem Q can be solved by iterative methods [2] . If in addition the Hopf-Lax 
operator A satisfies a causality property (see Lemma 1.3), then the fast marching algorithm [10] 
can be applied to decouple Q and solve it in a "single pass" , using a clever ordering of the set S7* . 

In order to make a sharp error analysis of this numerical scheme, we restrict as announced 
our attention to a very specific (and academic) setting. We consider the domain Q, = ]R d \ {0}, 
thus dVL = {0}, which will be discretized on the cartesian grid : f2* := 2Z d \{0}, <9f2* := {0}. We 
fix a matrix M G Sj" and we consider the constant riemannian metric defined by A4(z) = M 



(5) 



for all z G H . We denote by Dm the solution of p,h, which has the explicit expression 



\ z \\m, 



z G R d 



We use the variant introduced in [6] of the Hopf-Lax update operator A, which is based on 
a M-reduced mesh T (see Definition 1.1 in [B]), fixed in the rest of this paper. In other words 
T is a finite conforming mesh which satisfies the following properties 

(I) The union of the simplices T G T is a neighborhood of the origin. 

(II) The vertices of each simplex T G T lie on the lattice 2Z rf , and T has volume l/d\. 

(III) For each T G T, one of the vertices of T is the origin 0, and the others denoted by 
, Vd satisfy for all 1 < i,j < d the acuteness condition 



(vi,Vj) M > 0. 

For each map d : 7L d — > TR + U {oo} and each z G 7Z d , we define 

+ ^2 oti d(z + Vi) 



A(d, z) := min 

k,{oti),{vi) 



^ CiiVi 
Ki<k 



Ki<k 



(6) 



(7) 



where the minimum is taken among all 1 < k < d, all oti, ■ ■ ■ , a& G H+ such that aj+- • -+ak = 1, 
and all non-zero vertices v\, ■ ■ ■ ,Vk of a common simplex T G T ■ By convention, x oo = 0. 
For information, in the case of a non-constant riemannian metric Ai, a distinct A^(z)-reduced 
mesh T z is attached to each discrete point z and used to define the operator A(-, z), see [6]. 
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Our first result bounds the numerical error &m{z) — Dm(z), z G TL d , in terms of the geometry 
of the local mesh T ■ We introduce the bounding radius t"m(T) G K + , the radii r^(T) G H + 
and the angles 6m(T) G [0, tt] , for T G 7", defined by 

r&f(T) := max vmCT), tm{T) := max IUIIm, cos#j\/(T) := min ,, '^^5 — • (8) 
reT 2eT n,^er\{o} |M|m|M|m 

Note that the extrema defining rju(T) and 6m(T) are attained for vertices of T. We denote by 
fft+T the cone spanned by a simplex T G T : 

R + T := {rz; r G M+, * G T}. 

Theorem 1. If the Hopf-Lax update operator A is defined by ([?p, then the fixed point problem 

(9) 



d(*)=A(d,z), ze?Z d \{0}, 
d(0) = 

/ias a unique solution, which can be obtained "in one pass" using the fast marching algorithm. 
Denoting by <1m '■ TL d — > IR+ this solution, one has for all z G 7L d 

< d M (z) - D M (z) < dr M (T) (l + ln+ f^^y) ) , (10) 

where ln + (r) := maxjlnr, 0}. More precisely, if —z G iR+T for some T 6 T, i/ten 

< d A/ (z) - D M (z) < d (sin M (T)) 2 r M (T) (l + ln + (j^)) ■ ( n ) 



Inequality (10) states that the discrete output djvf of the fast marching algorithm overesti- 
mates the exact solution Dm by a logarithmic factor, which is at most proportional to bounding 
radius rjv/(T). This should not be a surprise, since the bounding radius rjv/(T) reflects the 
discretization step as seen by the matrix M, and since this algorithm relies on a first order 
discretization of the eikonal equation. The slightly sharper estimate (11) takes into account the 
angle width of the simplex T. 

Our second main result is an average estimate of the radius rjvf(T), when the M- reduced 
mesh T is constructed as described in [6], in dimension d < 3 and as described in Proposition 



3.2 in dimension 4. A key feature of this construction is that it guarantees a uniform bound on 



the mesh cardinality (#(T) < 6, 24 or 768 if d = 2, 3 or 4 respectively), and therefore a small 
complexity of the resulting numerical scheme, independently of the matrix M and thus of its 
anisotropy ratio k{M). 

We introduce the successive Minkowski minima [6] of a matrix M G S'j" , defined for 1 < i < d 

by 

Xi(M) := min{||ttj||Ar; (ui, • • • , ttj) G 7L d is free, and \\u\\\m < < |K||m}- (12) 



Note that \\{M) < ■ ■ ■ < A^(M) by construction. Choosing (ui, • • • , Ud) as the canonical basis, 
we obtain that \d(M) < ||M||a, an upper bound which is attained in the case of a diagonal 
matrix. We will use Minkowski's second theorem on successive minima ([5], p 199), a classical 
result of the geometry of numbers which states that for any M G 5j" 

2 -VdetM < Ai(M)---A d (M) < — VdetM, (13) 



d\u) d LU d 



4 



k = 30 



K= 100 



k = 1000 



3.0 
2.5 
2.0 
1.5 
1.0 
0.5 




0.0 



0.2 



0.4 



0.6 



o.x 




o.o 



0.2 



0.4 



0.6 



0.8 



IjJjjiiiJilliL 



0.0 0.2 0.4 0.6 0.8 



Figure 2: Graph of 9 \- > ^{RqDRq), where D is a diagonal matrix of entries (k, and 
Re £ O2 a rotation of angle 6 E [0, tt/4]. 



where cod denotes the volume of the d-dimensional euclidean unit ball. 

It follows from Corollary 1.8 in [6] that for any M G S^, 1 < d < 4, and any M-reduced 
mesh, one has 

Ad(M) < r M (T). (14) 



Combining this inequality with Minkowski's second theorem on successive minima (U3|, one 

1 1 — 

obtains the lower bound rd(7) > Xd{M) > Cd(det M)2i , where Cd = 2[d\oJd) 3 and denotes 

the volume of the d-dimensional unit ball. The next theorem shows that this lower estimate for 

rjVjf(T) is also an upper estimate, in an average sense, when the mesh T is constructed with our 

methods. 

Theorem 2. Let M E St. If 1 < d < 4 and T is a M-reduced mesh constructed as described 
in Proposition 1.9 or Proposition 1.10 of [6], or Proposition 



3. 2 of this paper, then 



r M (T) < K d \d(M). 



(15) 



with K 2 = 2, K 3 = 3, K 4 
of M and such that 



5. Furthermore for any d > 1 there exists a constant Cd independent 



X d {R T MR) dR < C d (det M)* 



(16) 



We denoted by Od the compact group of d x d orthogonal matrices, equipped with the canonical 
Haar probability measure. 

Let us immediately stress that the scaling factor (det M) 23 does not depend on the anisotropy 
of M : it is the same for the identity matrix and for a 2 x 2 matrix of eigenvalues A, A -1 . This 
factor only reflects the homogeneous scaling with M of the Minkowski minimum A^(M) (as well 
as the bounding radius rjw(T) and the exact or approximated distances Dm and djw). In view of 
this natural scaling, we may therefore limit our attention to matrices satisfying det(M) = 1, in 



which case the average upper estimate (16) is fully independent of M. In contrast the uniform 
upper bound ||-M||5 > A^(M) grows like a power of the anisotropy ratio, since 



K(M) d < 



\M\\ 



(det M) 



< k(M)~. 



As established in[2j and illustrated on Figure [2j such large values of A^(M) are statistically 
rare - in fact they correspond to pathological situations where an eigenvector associated to the 
small eigenvalue of M is equal or close to a small element of 7L d \ {0}, see jj2j Theorems [l] 
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and [2] therefore yield together a genuine result of non-linear approximation, which shows that 
the accuracy of the algorithm introduced in [6] does not degrade as the anisotropy ratio k{M) 
increases, at least in the case of a constant metric and in an average sense over all anisotropy 
orientations. 

We prove Theorems [T] and [2] respectively in |l] and ^2] A numerical experiment is presented 
in £j3j as well as the extension of the algorithm to dimension four. 



1 Error analysis of anisotropic fast marching 



This section is devoted to the proof of Theorem [TJ which is based on the notions of (discrete) 
sub-solution and super-solution of the fixed point discretization Q of the eikonal equation. We 
refer to [2] for more discussions on these notions. 
We allow discrete maps to take the value +00. 

Definition 1.1. A discrete map d : 7L d — > iR + U{oo} is a super-solution (resp. sub-solution) of 
the system if it satisfies 



d(z) > A(d,z), z£ffi d \{0}, 
d(0) = 0. 



yresp. 



d(*)<A(d,*), z£7Z d \{0}, 
d(0) = 0. 



(17) 



The exact continuous solution of the eikonal equation happens, in our specific setting, to be 
a discrete sub-solution, as shown in the next proposition. 

Proposition 1.2. The restriction of the exact solution D^/ (z) = \\z\\m, of the continuous eikonal 
equation to the lattice 7L d , is a discrete sub-solution of the system Q). 

Proof. Consider z G 7L d , 1 < k < d, non-negative coefficients (aj)f =1 , and vertices (vi)^ =1 of a 
common simplex T G T as they appear in the Hopf-Lax operator ([7]). Since a\ + • • • + = 1 
we obtain using first convexity and second the triangle inequality 



>M 



+ ^ a i v i 
Ki<k 



oii T) M {z + Vi) > Bj, 

l<i<k 

It follows that Dm(z) < A(Dm,z). Since Dm(0) = this concludes the proof. 



j > B M (z) - 






l<i<k 



M 



□ 



We next recall an argument often referred to as the Causality property, which is at the 
foundation of the fast marching algorithm. See [8] or [6j for a proof. 

Lemma 1.3 (Sethian Vladimirsky, 2000, [Sj ) . Let z G TL d and let d : 7L d ->■ M + U {00}. // 
A(d, z) < 00, then by compactness the minimum ([?]) defining the Hopf-Lax update A(d, z) is 
attained for some k G {1, ■ • • , d}, positive coefficients a\, • ■ ■ ,a,k > 0, a% + ■ • • + otk = 1, and 
non-zero vertices v\ , ■ ■ ■ ,Vk of a common simplex T G T . We then have for all 1 < i < k 



d(z + vi) < A(d,z). 



(18) 



The next proposition shows that any discrete super-solution d+ is larger than any sub- 
solution d_. A similar property is proved in [2], on a finite discrete domain, by studying the 
point where the difference (d+ — d_) reaches its minimum. Since the domain 7L d is infinite we 
cannot rely on this approach here, and we take advantage instead of the causality property. 
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Proposition 1.4. Let d+,d_ : 7L d — > M+ U {00} be respectively a discrete super- solution and 
sub-solution of the system |$|). Then d_ < d + on 7L d . 

Proof. We denote by r (resp. R) the maximum (resp. minimum) radius such that the ellipsoid 
{u G 1R d ; ||u||a/ < r} is contained in (resp. contains) the union of the elements of 7~, which is 
by assumption a compact neighborhood of the origin. Using first the fact that d+ is a super- 
solution, and second the definition ^ of A, we obtain that for each z G 7L d \ {0} there exists a 
vertex v of a simplex T G T such that 

d+(z) > A(d+, z)>r + d + (z + v). 



Since d+ takes non- negative values, and \\z + v\\m > \\ z \\m ~ R) it immediately follows that 
d+(z) > (r/i?) || z || . Hence the set {z G 7L d ] d + (z) < K} is finite for any constant K < 00, 
and it is therefore possible to sort the elements of {z G 7L d ; d + {z) < 00} into a sequence (z n )n>o 
ordered by increasing values : d + (z n ) < d + (z n+ i) for all n G 7L + . 

We prove that d_(z n ) < d+(z n ) by induction on n G 7L + . We first observe that d_(zo) = 
= d + (£o), since 20 = 0, and we next consider an arbitrary but fixed n G 7L + . In view of 
Lemma 1.3, there exists k G {1, • • • , d}, some positive coefficients (ctj)^ =1 summing up to 1, and 
ni, • • • , nk G 2Z+ such that 



A(d+,z n ) 



i<j<fc 



Ai 



+ 2^ "i d + , 

Ki<k 



Lemma 1.3 states in addition that d+(z ni ) < A(d+,.z n ) < d+{z n ), and therefore rij < n for all 
1 < i < k. It follows that d_(z n J < d + (z ni ) by induction, and therefore that A(d_,z n ) < 
A(d + , z n ). This implies d_(z n ) < d+(z n ), and concludes the proof of this proposition. □ 

Lemma 1.5. Let z G 7L d , and let T G T be such that —z G M+T. Denoting by vi, ■ ■ ■ ,Vd the 
non-zero vertices ofT, there exists non-negative integers /3\, ■ ■ ■ , (3d such that 

z + (3 1 vi + ■■■ + (3 d v d = 0. (19) 

Proof. Since the union of the elements of T is a neighborhood of the origin (property (I) of 
the mesh T), there exists as announced a simplex T such that — z G H+T. This inclusion 

implies the existence of non- negative reals such that z + P±v± H h fidX>d = 0, and we thus only 

have to show that /3±, ■ ■ ■ , (3d ar e integers. Since the vertices vi, ■ ■ ■ ,v d belong to 7L d , and since 
I det(^i, • • • , v d )\ = d\\T\ = 1 by assumption (property (II)), • • • , Vd) is a basis of the lattice 
7L d . The coefficients /3\, ■ ■ ■ , (3d are therefore integers, which concludes the proof. □ 

The variant studied here of the fast marching algorithm was shown in |6j to produce a 
solution of the system Q in a finite number of steps, in the case of a general continuous metric 
on a compact periodic domain. This algorithm involves local meshes T z for each discrete point 
z, which in the case of a constant metric are all equal to T ■ In order to adapt the proof of 
[6] to the case of a constant metric on an infinite domain, we briefly outline its main features. 
This algorithm produces a decreasing sequence (d n ) n >o of super-solutions of the system ^ : 
d n (^) < d m (z) for all n < m and all z G 7L d . The first element do of this sequence satisfies 
do(0) = and do(-z) = 00 for all z 7^ 0. Simultaneously, this algorithm produces an injective 
sequence of points (z n ) n >o which has the following properties : 

(i) d n (z n ) is the minimum value of d n on the set 7L d \ {z m ; m < n}. 

(ii) d n (z n ) = d m (z n ) = A(d m , z n ) for all m > n > 0. 

(iii) d n+ i takes finite values on the set {z G 7L d ; z n = z + v for some vertex v of T z }, and 
coincides with d ra outside of this set. 
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Proposition 1.6. The set Z := {z n ; n > 0} coincides with 7L d , and a solution d to the system 
^ is given by the decreasing limit, for all z G 7L d , 



d(z) : = lim d n (z). 



(20) 



Proof. We first observe that d n (z) > \\z\\m for all z G TL and all n > 0, since d n is a super- 
solution of the system ([9]), and zi-> a sub-solution. 

Let n > be arbitrary, and let v be an arbitrary vertex of T ■ Since d n+ i(z n — v)<oo (using 
property (hi) of the sequence (z n ) n >o, an d recalling that T = T z for a constant metric), the 
cardinality n' := E Z5 d ; ||x||m < d n+ i(z n — v)} is finite, which shows that z n — v = z m for 
some m < n' (using (i)). 

As a result for any z G Z and any vertex v of T we have z — v E Z. Since the minimum 
value of <io is attained at the origin only, we have zq = and therefore G Z. Since any z G 7L d 
can be written under the form (19), it follows that Z = 7L d . In view of (ii), the discrete map 
defined by ( 20 ) is a solution of the system Q , which concludes the proof. □ 



The fast marching algorithm is regarded as a "one pass algorithm" because the sequence 
(z n )n>o is injective. In the case of a constant metric this implies in particular, using property 
(iii), that for any z G 7L d , the cardinality of {n > 0; d n+ \(z) / d n (z)} is bounded by the number 
of vertices of T ■ Proposition 1.6 establishes the first part of Theorem [TJ while Propositions 1.2 
and 1.4 imply the lower estimates for d*/ in (10) and (11). 



1.1 Construction of a discrete super-solution 



We construct in this section an explicit discrete super-solution d + of the system ([9]). Using 
its expression, and the fact that super-solutions are larger than solutions, we obtain the upper 
estimates for announced in Theorem [TJ ( |10[ ) and (11), which concludes its proof. 

We consider an arbitrary simplex T G T, which is fixed throughout this section, and we 
define a discrete map d + : 7L d — > IR+ U {oo} as follows. If z G 7L d is such that —z ^ 1R+T, then 
we set d + (z) := oo. Otherwise we write z under the form (19), and we define 



d+(z) := \\z\\ M + { S m9 M {T)) 2 £ s(A)|| Vi\\M, where s((3) := 



Ki<d 



E 

l<fc</9 



The summation appearing in d+ can be bounded as follows : 



s(f3i)\\vi\\ M < ^2 iMUffi+in 



Ki<d 



Ki<d 



\\z\\m 

\Vi\\M 



< dr M {T) 1+ln 



\\Z\\M 

r M {T) 



(21) 



(22) 



where we used in the first inequality the upper bound < \\z\\M/\\vi\\M (due to the acuteness 
condition Q), and Ylk=i ^ l — 1 + m+ (/3)- For the second inequality we used the growth of 
r i— > rln(l + ln + (A/r)) on 1R^_, for any fixed A > 0. Let us assume for a moment that d+ is a 
super-solution, which implies that d+ > djv/. Then combining (22) and (21) we obtain the upper 



estimates in (11). Since the simplex T G T is arbitrary, we also obtain (10), which concludes 



the proof of Theorem [T] 

In order to prove that d+ is a super-solution, we need a preliminary lemma, on the Taylor 
development of the || • ||m norm, under an acuteness condition. 
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Lemma 1.7. Let z,v G R d \ {0} and 6*o G [0, tt/2] be such that 



(v, z - v) M > cos(0 o )|M|m||2 - v\\ M - 



(23) 



Then 



\z - v\\ M < \\z\\m 



{z,v) m , (sin6>o 



2 M 



+ 



\M 



\ Z \\M 



Proof. Since the vectors v and z—v form an angle smaller than 9$ ( 23 ) , this is also the case for the 



vectors v and z— tv = z—v+(l—t)v, for £ £ [0, 1]. Therefore («, z—tv) m > cos(0o)IM|m||-z — £^||mj 
for £ G [0, 1]. On the other hand we have ||.z||m|M|a/ > (#) > (z — v, v)m + IMIm — IMIm' 
hence [|z[|m > H^Ha^ which implies that ||« — tv\\M > IH|a/ — i[| v ||m > (1 — £)||z||m, for £ G [0, 1]. 
The Taylor formula with integral rest, applied to the smooth function £ i-> ||z — £v||a/, yields 



ili 



+ 



ili 



(z — tv, v) 



M 



tV\\M 



(1 - t)dt. 



Injecting the lower bounds (z — tv, v) m > cos(#o)IM|a/ \\z — £u||iw, and \\z — £u||m > (1 — £) IHImj 
we obtain that the integrated term is bounded above by (sin ) 2 1 1 ^ 1 1 / 



z\\m- This concludes 

the proof of this lemma. □ 

Proposition 1.8. The map d+ defined by (21) is a discrete super- solution of Q). 

Proof. As required, we have d+(0) = 0. We thus consider an arbitrary z G 7L d \ {0}, and we 
show in the following that d+(z) > A(d+, z). There is nothing to prove if d+(z) = +oo, and we 

0, where fa G 7L + for 1 < i < d 



may therefore write z under the form ( 19 ) : z + Yli=i fii v i 
and v\, - ■ ■ ,Vd are the vertices of T. 



For all 1 < i < d such that fa > 0, we apply Lemma 1.7 to z, v := —V{, and the angle &m(T). 



We obtain, since \\z\\m > fa\\vi\\M (due to the acuteness condition (|6 

{-Z,Vi) M 



z M 



\Z + Vi\\M > 



z M 



(sin 8 M (T)Y 



^> ^ -(smMT)f IKiM 



PA/ 



Hence d+(z) — d+(z + Vj) > {—z,Vi)M/\\z\\M, for all 1 < 2 < d such that fa > 0. Defining 
ai := fa/ (fa H h ^d), in such way that Xy«=i ^i^i is positively proportional to — z, we obtain 



E 

Ki<rf 



+ ^2 a i d+(z + Vi) < 



M 



Ki<d 



E 

Ki<d 



CtiVo 



Ai 



a<i<k 



z M 



d + (z). 



ili 



This establishes that d + (z) > A(d+,z), which concludes the proof. 



□ 



2 Average value of the outer radius tm{T) 

This section is devoted to the proof of Theorem [2j which begins with the upper bound (15) on 
r Ai(T). An inspection of the construction of the mesh T in Proposition 1.9 or Proposition 1.10 
of |6j , in dimension d = 2 or 3 respectively, shows that its vertices have the form respectively 

£iui + e%ui, or s\u\ + e 2 U2 + £3^3, 

where || JWf = K(M) and E{ G {—1,0,1}, for all 1 < % < d. The upper bound on tm(T) 
thus immediately follows from the triangle inequality and Ai(M) < • • • < A^(M). In the four 
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dimensional case, which is covered by Proposition 3.2 of this paper, the vertices also have the 
form e\U\ + £2^2 + £3^3 + E4U4, where three of the \si\ are bounded by 1 and the remaining one 
is bounded by 2. The announced inequality (15) again follows from the triangle inequality. 



We next turn to the proof of (16), in which the dimension d > 1 is arbitrary. The orthogonal 
group Od is equipped with the canonical Haar probability measure. We denote by R a random 
orthogonal matrix, and by V the probability of an event in Od- 

For each matrix M G S^, we denote by < ^i(M) < • • • < Ud(M) the square roots of the 
successive eigenvalues of M. In particular v\{M) = || A/ -1 1| — 2 ? and Ud{M) = ||M||2. We denote 
by \u\ := V u T u the euclidean norm of a vector u G H d . 

Lemma 2.1. Let M G Sf be such that det(M) = 1. Then for any u G M d such that \u\ = 1 
one has 

V(\\Ru\\ M <S)< C(d)5 d ~ 1 iy 1 {M), (24) 
where C(d) = 2 d /uid, and u>d is the volume of the d- dimensional euclidean unit ball. 
Proof. We denote by B := {x G lR d ; \x\ < 1} the euclidean unit ball, and we introduce the set 

B := {x G B; \\x\\ M < 5\x\}. 

Since the Haar measure on Od is invariant under the action of rotations, its image by R >— > Ru is 
the uniform probability on the euclidean sphere. The probability estimated in (24) is therefore 
equal to |£?o|/|-B|. The invariance of the Haar measure under rotations also allows us to assume 
without loss of generality that M is a diagonal matrix. For notational simplicity we denote 
Vi := fj(M), for 1 < i < d, and we observe that v\ ■ ■ ■ v& = V det M = 1. 
If x = (xi, • • • ,Xd) G Bo, then 



22, 1 22^r-2ii2^c-2 

v{x\ + • • • + Vd x d — ° \ x \ ^ v ) 

hence 

x G [-1, 1] x [S/u 2 , 5/v 2 ] x • • • x [-8/v d , 5/u d \. 
It follows that \Bq\ < 2 d 5 d ~ 1 / (v 2 ■ ■ ■ ^d) = 2 d 5 d ~ l v\, which concludes the proof. □ 
Our next result is an estimate in probability of the first Minkowski minimum of M : 

Ai(M) := min 

u£2Z d \{0} 

Corollary 2.2. Let M G Sj be such that det(M) = 1. Then for each 5 > 

V( Ai(i? T Mi?) < 5 ) < C'(d) 5 d , (25) 
where C'(d) := 2d3 d ~ 1 C(d), and C{d) denotes the constant from Lemma 



2.1 



Proof. We introduce the collection of points E := {u G 7L d ; < |n| < 5/ui(M)}, and we observe 
that \\Ru\\ M > $ for any u G 7L d \(EU {0}) and any R G O d . We obtain 

V{X 1 (R T MR) <5) < V(\\Ru\\ M < 5) 

u£ZL d \{0} 

= Y,V(\\Ru\\ M <5) 

< C{d)u 1 {M)Y J {5/\u\) d -\ (26) 
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where in the third line we applied Lemma 2.1 to the normalized vectors u G E. 

In order to estimate ^ ue g((5 ^/|u|) a! ~ 1 , we introduce the || • ||oo norm on JR d defined by 
||(ui,-- - ,w,i)||oo := max{|ui|, ■ ■ ■ , |ud|}. We observe that ||w||oo < and that for each positive 
integer k there exists exactly (2k + l) d — (2k — l) d elements u G 7L d which satisfy ||w||oo = k. 
Therefore 

(2k + l) d - (2k - l) d 



E 



il-d 



E 

Q<k<8/v 1 {M) 

Remarking that (2k + l) d -(2k- l) d < d(2k + if- 1 ((2k + 1 
bound the right hand side of (27) by (5 j 1 v\(M))2d?> d ~ 1 . Combining this estimate with (26), we 
conclude the proof. 



k d-i 
)-(2k 



(27) 



1)) < 2d(3k) d - 1 , we 



□ 



The right part of Minkowski's second theorem on successive minima (13) implies 

Ad(M) < 

For any matrix M G Sj" such that det(M) 
X d (R T MR)dR 



2 d 



< 



o,, 



2 d VdetM 
^Ai(M) rf -!- 

1, we thus obtain 

dR 



(d 



Xx^MRy- 1 

roc 

- 1) / V( Ai (R T MR) < 5 ] 
Jo 



dS 



< (d-1) 



50 d6 



(d- l)C'(d) + 1. 



This concludes the proof of (16) in the case of a matrix M of determinant 1. 

The general case of an arbitrary determinant can be reduced to the case det(M) = 1 by 
homogeneity. Indeed the determinant is homogeneous : det(> 2 M) = fi 2d det(M), for any \i > 0; 
and so are Minkowski's minima : Xi(fj, 2 M) = fj,Xi(M) (this is an immediate consequence of the 
homogeneity of the norm : 



u 



p. M = ^||u||m for any «£]R 



3 Numerical experiments, and extension to dimension 4 

We discuss in the first part of this section some two-dimensional numerical experiments, which 
illustrate our two main results Theorem [T] and Theorem [2] These numerical experiments are con- 
ducted in the restrictive setting which is the framework of these results: a constant riemannian 
metric. The reader interested in benchmarks closer to applications is referred to [6]. The second 
part of this section is devoted to the extension of our variant of the fast marching algorithm to 
four dimensional domains. 



We conducted some numerical experiments, in which the anisotropic eikonal equation as- 
sociated to a constant riemannian metric M = M G is discretized on the finite grid 
SI* := {-500,- •• ,500} 2 \ {(0,0)} and solved using three different methods. The matrix M 
has eigenvalues 1/10 and 10, thus anisotropy ratio k(M) = 10, and the eigenvector associated to 
the eigenvalue 1/10 is (1,0.6). The fast marching algorithm (FM) is applied with two different 
M-reduced meshes 71 and T2 presented on Figure [3] (top and bottom respectively) , of cardinality 
6 and 10 respectively. The mesh 71 is constructed as described in Proposition 1.9 in [6]. The 
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Figure 3: The two meshes 7i (top) and T2 (bottom). Portion of the domain where our error 
analysis of the fast marching algorithm applies (thick, dashed, red arrows) for T\ (center) and 
T2 (right). 



mesh 7i is obtained using a different strategy, which will be the object of a future article : 
an iterative subdivision procedure which stops when the mesh satisfies the acuteness condition 
Q. The results are compared with with the Iterative Iterative Gauss Siedel (IGS) algorithm 
introduced in [2J. 





IGS 


FM 71 


fm r 2 


|| D 


- d L°°(n») 


25.2 


2.69 


0.803 


D 


-d||l,°o(Ql) 


25.2 


1.25 


0.803 




Timing 


2.0s 


3.6s 


13.1s 



The computation times were obtained on a 2.4Ghz core 2 duo laptop. Surprisingly, the IGS is 
the fastest algorithm on this test case, but it is also the most imprecise. Indeed the numerical 
error (on the first line of the table) is reduced by a factor 9 when one replaces the IGS with 
the fast marching algorithm with the mesh 71, and by an additional factor 3 with the mesh 72- 
For comparison, the original isotropic fast marching algorithm, with M' = Id and on the same 
domain, yields a numerical error of 2.1, which confirms that our variant of the fast marching 
algorithm did not suffer from the anisotropy of M , contrary to the IGS. 

The numerical error is defined as the maximum absolute value of the difference D — d between 
the exact solution of the eikonal equation D(z) = ||z||mj an d a discrete approximation d produced 
by the algorithm of interest. The maximum is taken either on the whole finite grid fi* for the 
first line of the table, or for the second line on a subset Ql which is described in the next remark. 

Remark 3.1. The error analysis presented in Theorem [i] was obtained on the infinite grid 
K 2 \{(0,0)} ; instead of the finite gridtl* = {— 500, ••• , 500} 2 \{(0, 0)} ; and is not automatically 
valid on the whole $7* for the following reason. Consider two non-zero vertices u,v of a triangle 
T in the M-reduced mesh T , and a, (3 G 7L + such that z := —{au + f3v) belongs to the discrete 
domain. We used in Proposition 1.8 that — (a'u + (3'v) also belongs to the discrete domain, for 
all0<a' <a and all < /?' < 0. 

This property does not hold for the mesh T\ and the finite grid f2*, but does hold for the mesh 
7~2 on the same grid, see Figure^ (center and right respectively). More precisely, for the meshTi, 
this property holds on a subset 0* which contains approximately 37% of the elements ofQ*, and 
is illustrated on Figure [3] (center; thick, dashed, red arrows). As expected, and illustrated on the 
second line of the table, the numerical error of the fast marching algorithm with the M-reduced 
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mesh 71 is significantly reduced on this subset, while it is unchanged with the mesh T2 or the IGS. 



This rest of this section is devoted to the proof of the following proposition, which extends 
the variant of the fast marching algorithm presented in [B] to four dimensional domains. 

A M-reduced basis, where M E S£, 1 < d < 4, is a collection (u±, • • • , U4) of elements of 7L d 
satisfying | det(iti, • • • ,Ud)\ = 1 and ||itj||M = Aj(M). Such a basis exists and can be obtained 
with an algorithm of complexity 0(ln k(M)), if one regards the elementary operations among 
reals (+, — , x,/) as of unit complexity, see (HE]. 

Proposition 3.2. Let M 6 and let (u%, 112, u 3 , 114) be a M-reduced basis. Consider the 
collection T of simplices which have the following vertices 

(0, vx, V1+V2, V1 + V2 + V3, 2v\ + V2 + v 3 + V4), (28) 

(0, V\ + V2, V1+V2 + V3, Vl + V2 + V 3 + 1> 4 , 1v\ +^2 + ^3+ V4), (29) 

where (i>i, i>2, i>3, i> 4 ) runs over all permutations of (emi, E2U2, £31x3, £41*4), with arbitrary signs 
(e%, 82, £3, £4) € { — 1, 1} 4 - Then T is a M-reduced mesh, of cardinality 768. 

The cardinality of #(7~) is indeed 2 x 4! x 2 4 = 768, where the factor 2 stands for the two 
types of simplices, 4! is the number of permutations of a four elements set, and 2 4 is the number 
of choices for the signs (ex, £2, £3, £4)- The fact that T is a conforming mesh, and that the 
union of its elements is a neighborhood of the origin, is easily checked when (u±, U2, U3, U4) is 
the canonical basis of TR d . It thus holds for an arbitrary basis (u\, 112, 113, U4) of ]R rf by change 
of variables. Observing that 

|det(ui, V1+V2, vi + V2 + v 3 , 2v\ + v 2 + v 3 + v 4 )| = | det(«i, u 2 , u 3 , Ui)\, 
I det(fi + v 2 , V1 + V2 + v 3 , ^1+^2 + ^3 + t> 4 , 2v\ +V2 + V3 + V4)\ = I det(«i, u 2 , u 3 , 1x4) | , 

and recalling that | det(ni, 112, u 3 , Ui)\ = 1, by definition, we obtain that the volume of each 
element of T is 1/4!. In order to conclude the proof, it only remains to establish the acuteness 
property 

We now recall a property of the M-reduced basis (u%,U2, 143,144), established in Proposition 
1.6 of [BJ. For any integer combination z of the elements of the basis distinct from Uj, in other 
words z = a.\U\ + • • • + OLi-\Ui-\ + cti + \Ui + i + • • • + «4U 4 , where a%, ■ ■ ■ , a^-i, on+i, • • • , a 4 ETL, 
one has 

2\(z, Ui )M\<\\z\\ 2 M . (30) 

The following lines establish that the angle (it, v) m formed by any pair of non-zero vertices 
u, v of any simplex T G T is non- negative. For that purpose we distribute by bi-linearity the 
terms of the scalar product (it, v)m, and we group them within square brackets. It easily follows 



from (30) that the sum of the contents of each square bracket is non-negative. In the following 



three cases, this inequality should be applied with z = v\. 



v\ + v 2 )m = [ IK Hat + {vi,v 2 )m]- 
{vi, vx + V2 + v 3 ) M = [IKIIm + {vi,v 2 )m + {vi,v 3 ) M \- 
{vi,2v l + V2 + v 3 + v i ) M = [2||«i||m + {vi,v 2 )m + {vi,v 3 )m + {vx^v^m. 
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In the next three cases, one should apply (30) with z = V\, z = v\ + v 2 or z = v\ + v 2 + v 3 . 



(vi + v 2 , vi+v 2 + v 3 ) M 

(Vl + V 2 , 2vi + V 2 + V 3 + V4) M 



(Vi + V 2 + V 3 , 2l>l + V 2 +V 3 + V4) 



M 



[ \\vi + V 2 \\m + {vi + V 2 ,V 3 ) M ]- 

[ IK + v 2\\ll + ( v l + V 2 ,V 3 )m + (Vl + V 2 ,V A ) M ' 
+ [ \\ v i\\m + («1,«2)m]- 

[ Ibl + V 2 + V 3 \\ 2 M + (Vl + V 2 + V 3 ,V4)m] 



+ [ IMIm + (vx,v 2 )m + {vi,v 3 ) M ]- 

The non-negativity of the above expressions shows as announced that any two non-zero vertices 
u,v of any simplex of the form (28) form a non-negative scalar product. We obtain a similar 
conclusion for simplices of the form (29) by checking that the following scalar products are non- 
negative. One should apply (30) with z = v±, z = v% + v 2 , z = v\ + v 2 + v 3 or z = v 2 + 1*3 + 1*4. 



(Vl + V 2 , Vl + V 2 + V 3 + V 4 ) M 
{Vl +V2+ V 3 , Vl + V 2 + V 3 + V A ) M 
(2t>i + V2 + V 3 + U4, Vl + U 2 + «3 + V4) M 



[ Ibl + V 2 \\m + («1 + V 2 , V 3 ) M + {Vl + V 2 ,V4) M . 

[ \\V1 +V 2 + V 3 \\ 2 M + (V! + V 2 + V 3 ,V^) M ]- 

[ \\v 2 +V 3 + V4\\m + 2{v 2 + V 3 + Vi, Vi) ] 

+ [2\\vi\\ 2 M + {vi,v 2 )m + {vi,v 3 ) M + (V1,V4) M , 



Conclusion 

We have continued in this paper the analysis of a variant of the fast marching algorithm, intro- 
duced in [6], which is particularly efficient in the context of large anisotropies. The computational 
complexity of this method was known to be largely insensitive to anisotropy. This paper estab- 
lishes, in the special case of a constant metric on the domain H rf \ {0}, 1 < d < 4, and in an 
average sense over all anisotropy orientations, that the numerical accuracy not affected either 
by anisotropy, however pronounced. A two-dimensional numerical experiment confirms this ex- 
cellent accuracy, and an extension of the algorithm to four dimensional domains is presented. 
Future work will be devoted to the error analysis in the case of a general continuous metric, and 
to the application of this algorithm to medical image and data analysis. 
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